clear all
set more off

use plos_data_individual.dta

* Figure 2 *

* Age
qui: sum age4 if witchcraft!=.
local N: di %7.0fc r(N)
qui: tab isocode if age4!=. & witchcraft!=.
qui: local C = r(r)
graph bar witchcraft, over(age4, label(labsize(large) labgap(small))) ///
note("Sample of `N' persons from `C' countries", pos(12) size(large)) ///
ytitle("Share of witchcraft believers") ylabel(0(0.1)0.5, labsize(large)) ysc(r(0 0.55)) scheme(gap_axes) graphregion(color(white) margin(0 0 0 1)) bgcolor(white) legend(off)
graph export "Figures\Fig2_1.pdf", replace

* Gender
qui: sum gender if witchcraft!=.
local N: di %7.0fc r(N)
qui: tab isocode if gender!=. & witchcraft!=.
qui: local C = r(r)
graph bar witchcraft, over(gender, label(labsize(large) labgap(small))) ///
note("Sample of `N' persons from `C' countries", pos(12) size(large)) ///
ytitle("Share of witchcraft believers") ylabel(0(0.1)0.5, labsize(large)) ysc(r(0 0.55)) scheme(gap_axes) graphregion(color(white) margin(0 0 0 1)) bgcolor(white) legend(off)
graph export "Figures\Fig2_2.pdf", replace

* Location of residence
qui: sum urban if witchcraft!=.
local N: di %7.0fc r(N)
qui: tab isocode if urban!=. & witchcraft!=.
qui: local C = r(r)
graph bar witchcraft, over(urban, label(labsize(large) labgap(small))) ///
note("Sample of `N' persons from `C' countries", pos(12) size(large)) ///
ytitle("Share of witchcraft believers") ylabel(0(0.1)0.5, labsize(large)) ysc(r(0 0.55)) scheme(gap_axes) graphregion(color(white) margin(0 0 0 1)) bgcolor(white) legend(off)
graph export "Figures\Fig2_3.pdf", replace

* Education
qui: tab educ if witchcraft!=.
local N: di %7.0fc r(N)
qui: tab isocode if educ!=. & witchcraft!=.
qui: local C = r(r)
graph bar witchcraft, over(educ, label(labsize(large) labgap(small))) ///
note("Sample of `N' persons from `C' countries", pos(12) size(large)) ///
ytitle("Share of witchcraft believers") ylabel(0(0.1)0.5, labsize(large)) ysc(r(0 0.55)) scheme(gap_axes) graphregion(color(white) margin(0 0 0 1)) bgcolor(white) legend(off)
graph export "Figures\Fig2_4.pdf", replace

* Personal economic situation
qui: sum econper if witchcraft!=.
local N: di %7.0fc r(N)
qui: tab isocode if econper!=. & witchcraft!=.
qui: local C = r(r)
graph bar witchcraft, over(econper, label(labsize(large) labgap(small))) ///
note("Sample of `N' persons from `C' countries", pos(12) size(large)) ///
ytitle("Share of witchcraft believers") ylabel(0(0.1)0.5, labsize(large)) ysc(r(0 0.55)) scheme(gap_axes) graphregion(color(white) margin(0 0 0 1)) bgcolor(white) legend(off)
graph export "Figures\Fig2_5.pdf", replace

* Household size
qui: sum hhsize3 if witchcraft!=.
local N: di %7.0fc r(N)
qui: tab isocode if hhsize3!=. & witchcraft!=.
qui: local C = r(r)
graph bar witchcraft, over(hhsize3, label(labsize(large) labgap(small))) ///
note("Sample of `N' persons from `C' countries", pos(12) size(large)) ///
ytitle("Share of witchcraft believers") ylabel(0(0.1)0.5, labsize(large)) ysc(r(0 0.55)) scheme(gap_axes) graphregion(color(white) margin(0 0 0 1)) bgcolor(white) legend(off)
graph export "Figures\Fig2_6.pdf", replace

* Religious affiliation
qui: tab religion3 if witchcraft!=.
local N: di %7.0fc r(N)
qui: tab isocode if religion3!=. & witchcraft!=.
qui: local C = r(r)
graph bar witchcraft, over(religion3, label(labsize(large) labgap(small))) ///
note("Sample of `N' persons from `C' countries", pos(12) size(large)) ///
ytitle("Share of witchcraft believers") ylabel(0(0.1)0.5, labsize(large)) ysc(r(0 0.55)) scheme(gap_axes) graphregion(color(white) margin(0 0 0 1)) bgcolor(white) legend(off)
graph export "Figures\Fig2_7.pdf", replace

* Religiosity
qui: sum religiosity if witchcraft!=.
local N: di %7.0fc r(N)
qui: tab isocode if religiosity!=. & witchcraft!=.
qui: local C = r(r)
graph bar witchcraft, over(religiosity, label(labsize(large) labgap(small))) ///
note("Sample of `N' persons from `C' countries", pos(12) size(large)) ///
ytitle("Share of witchcraft believers") ylabel(0(0.1)0.5, labsize(large)) ysc(r(0 0.55)) scheme(gap_axes) graphregion(color(white) margin(0 0 0 1)) bgcolor(white) legend(off)
graph export "Figures\Fig2_8.pdf", replace

* Belief in god
qui: sum god if witchcraft!=.
local N: di %7.0fc r(N)
qui: tab isocode if god!=. & witchcraft!=.
qui: local C = r(r)
graph bar witchcraft, over(god, label(labsize(large) labgap(small))) ///
note("Sample of `N' persons from `C' countries", pos(12) size(large)) ///
ytitle("Share of witchcraft believers") ylabel(0(0.1)0.5, labsize(large)) ysc(r(0 0.55)) scheme(gap_axes) graphregion(color(white) margin(0 0 0 1)) bgcolor(white) legend(off)
graph export "Figures\Fig2_9.pdf", replace


* Table 1 *

eststo T1_1: xi: dprobit witchcraft age gender i.isocode, clu(isocode)
eststo T1_2: xi: dprobit witchcraft age gender i.educ i.isocode, clu(isocode)
eststo T1_3: xi: dprobit witchcraft age gender i.educ i.econper i.isocode, clu(isocode)
eststo T1_4: xi: dprobit witchcraft age gender i.educ i.econper i.hhsize3 i.isocode, clu(isocode)
eststo T1_5: xi: dprobit witchcraft age gender i.educ i.econper i.hhsize3 urban i.isocode, clu(isocode)
eststo T1_6: xi: dprobit witchcraft age gender i.educ i.religion3 i.religiosity i.isocode, clu(isocode)
eststo T1_7: xi: dprobit witchcraft age gender i.educ i.religion3 god i.isocode, clu(isocode)
eststo T1_8: xi: dprobit witchcraft age gender i.educ i.econper i.hhsize3 urban i.religion3 i.religiosity i.isocode, clu(isocode)

estadd scalar cnt_cls=e(N_clust): *

esttab T1* using "Tables\Table1.tex", drop(_Iisocode* _cons) label replace b(3) se(3) booktabs alignment(D{.}{.}{3.5}) ///
se star(^{*} 0.1 ^{**} 0.05 ^{***} 0.01) margin title() nomtitles nonotes nodepvar ///
stats(N cnt_cls, fmt(%9.0f %9.0f) labels("Observations" "Countries"))

eststo clear

* Table A.2 of S1 Appendix *

tab educ, gen(educ_s)
tab econper, gen(econper_s)
tab hhsize3, gen(hhsize3_s)
tab religion3, gen(religion3_s)
tab religiosity, gen(religiosity_s)

gen age10=age*10
estpost tabstat witchcraft age10 gender urban god educ_s* econper_s* hhsize3_s* religion3_s* religiosity_s* if witchcraft!=., c(stat) stat(mean sd min max n)

esttab using "Tables\TableA2.tex", replace cells("mean(fmt(%-9.3gc)) sd(fmt(%-9.3gc)) min max count") nonumber ///
  nomtitle nonote noobs label collabels("Mean" "St. dev." "Min" "Max" "Obs.") booktabs

drop *_s* age10
eststo clear

* Table B.1 of S1 Appendix *

encode isocode, gen(isocode_num)

eststo TB1_1: reg witchcraft age gender i.isocode_num, clu(isocode_num)
eststo TB1_2: reg witchcraft age gender i.educ i.isocode_num, clu(isocode_num)
eststo TB1_3: reg witchcraft age gender i.educ i.econper i.isocode_num, clu(isocode_num)
eststo TB1_4: reg witchcraft age gender i.educ i.econper i.hhsize3 i.isocode_num, clu(isocode_num)
eststo TB1_5: reg witchcraft age gender i.educ i.econper i.hhsize3 urban i.isocode_num, clu(isocode_num)
eststo TB1_6: reg witchcraft age gender i.educ i.religion3 i.religiosity i.isocode_num, clu(isocode_num)
eststo TB1_7: reg witchcraft age gender i.educ i.religion3 god i.isocode_num, clu(isocode_num)
eststo TB1_8: reg witchcraft age gender i.educ i.econper i.hhsize3 urban i.religion3 i.religiosity i.isocode_num, clu(isocode_num)

estadd scalar cnt_cls_lpm=e(N_clust): *

esttab TB1* using "Tables\TableB1.tex", drop(*.isocode_num 1.* _cons) label replace b(3) se(3) booktabs alignment(D{.}{.}{3.5}) ///
se star(^{*} 0.1 ^{**} 0.05 ^{***} 0.01) title() nomtitles nonotes nodepvar ///
stats(N cnt_cls_lpm, fmt(%9.0f %9.0f) labels("Observations" "Countries"))

drop isocode_num
eststo clear

* Table B.2 of S1 Appendix *

eststo TB2_1: xi: dprobit witchcraft age gender i.wave, clu(isocode)
eststo TB2_2: xi: dprobit witchcraft age gender i.educ i.wave, clu(isocode)
eststo TB2_3: xi: dprobit witchcraft age gender i.educ i.econper i.wave, clu(isocode)
eststo TB2_4: xi: dprobit witchcraft age gender i.educ i.econper i.hhsize3 i.wave, clu(isocode)
eststo TB2_5: xi: dprobit witchcraft age gender i.educ i.econper i.hhsize3 urban i.wave, clu(isocode)
eststo TB2_6: xi: dprobit witchcraft age gender i.educ i.religion3 i.religiosity i.wave, clu(isocode)
eststo TB2_7: xi: dprobit witchcraft age gender i.educ i.religion3 god i.wave, clu(isocode)
eststo TB2_8: xi: dprobit witchcraft age gender i.educ i.econper i.hhsize3 urban i.religion3 i.religiosity i.wave, clu(isocode)

estadd scalar cnt_cls_wav=e(N_clust): *

esttab TB2* using "Tables\TableB2.tex", drop(_Iwave* _cons) label replace b(3) se(3) booktabs alignment(D{.}{.}{3.5}) ///
se star(^{*} 0.1 ^{**} 0.05 ^{***} 0.01) margin title() nomtitles nonotes nodepvar ///
stats(N cnt_cls_wav, fmt(%9.0f %9.0f) labels("Observations" "Countries"))

eststo clear
